import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mplt
import seaborn as sns
import arviz as az
import datetime

import fit_arima

import open_data
df_temp = open_data.open()
df = df_temp[0]

ritorno = fit_arima.compute(2, 1, catene=4, samples_per_chain=2500, burnin=1000)

19:49:51 - cmdstanpy - INFO - CmdStan start processing

chain 1 |          | 00:00 Status

chain 2 |          | 00:00 Status

chain 3 |          | 00:00 Status

chain 4 |          | 00:00 Status

20:05:56 - cmdstanpy - INFO - CmdStan done processing.
20:05:56 - cmdstanpy - WARNING - Non-fatal error during sampling:
Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/home/br1/PythonProjects/PythonStats/git_bayesian_copie/ARIMA/OnlyCovariates/code.stan', line 157, column 4 to column 34)
	Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/home/br1/PythonProjects/PythonStats/git_bayesian_copie/ARIMA/OnlyCovariates/code.stan', line 157, column 4 to column 34)
	Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/home/br1/PythonProjects/PythonStats/git_bayesian_copie/ARIMA/OnlyCovariates/code.stan', line 157, column 4 to column 34)
	Exception: code_model_namespace::log_prob: phi[2][1] is -nan, but must be greater than or equal to -1.000000 (in '/home/br1/PythonProjects/PythonStats/git_bayesian_copie/ARIMA/OnlyCovariates/code.stan', line 87, column 2 to column 47)
	Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/home/br1/PythonProjects/PythonStats/git_bayesian_copie/ARIMA/OnlyCovariates/code.stan', line 157, column 4 to column 34)
Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/home/br1/PythonProjects/PythonStats/git_bayesian_copie/ARIMA/OnlyCovariates/code.stan', line 157, column 4 to column 34)
	Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/home/br1/PythonProjects/PythonStats/git_bayesian_copie/ARIMA/OnlyCovariates/code.stan', line 157, column 4 to column 34)
	Exception: code_model_namespace::log_prob: theta[1] is -nan, but must be greater than or equal to -1.000000 (in '/home/br1/PythonProjects/PythonStats/git_bayesian_copie/ARIMA/OnlyCovariates/code.stan', line 86, column 2 to column 36)
	Exception: code_model_namespace::log_prob: theta[1] is -nan, but must be greater than or equal to -1.000000 (in '/home/br1/PythonProjects/PythonStats/git_bayesian_copie/ARIMA/OnlyCovariates/code.stan', line 86, column 2 to column 36)
Exception: code_model_namespace::log_prob: phi[1][2] is -nan, but must be greater than or equal to -1.000000 (in '/home/br1/PythonProjects/PythonStats/git_bayesian_copie/ARIMA/OnlyCovariates/code.stan', line 87, column 2 to column 47)
	Exception: code_model_namespace::log_prob: phi[1][5] is 1, but must be less than or equal to 1.000000 (in '/home/br1/PythonProjects/PythonStats/git_bayesian_copie/ARIMA/OnlyCovariates/code.stan', line 87, column 2 to column 47)
Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/home/br1/PythonProjects/PythonStats/git_bayesian_copie/ARIMA/OnlyCovariates/code.stan', line 157, column 4 to column 34)
	Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/home/br1/PythonProjects/PythonStats/git_bayesian_copie/ARIMA/OnlyCovariates/code.stan', line 157, column 4 to column 34)
	Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/home/br1/PythonProjects/PythonStats/git_bayesian_copie/ARIMA/OnlyCovariates/code.stan', line 157, column 4 to column 34)
	Exception: code_model_namespace::log_prob: phi[1][1] is -nan, but must be greater than or equal to -1.000000 (in '/home/br1/PythonProjects/PythonStats/git_bayesian_copie/ARIMA/OnlyCovariates/code.stan', line 87, column 2 to column 47)
Consider re-running with show_console=True if the above output is unclear!

/home/br1/PythonProjects/PythonStats/ARIMA_last_computation/OnlyCovariates/fit_arima.py:41: FutureWarning: Indexing with a float is deprecated, and will raise an IndexError in pandas 2.0. You can manually convert to an integer key instead.
  y_missing_list.append({'Stazione':df.columns[v[i]],'Data':df.index[u[i]],'Samples':inference_data.posterior.w.values[:,:,i].reshape(catene*samples_per_chain)})

Tempo di computazione: 16:05

res = az.loo(ritorno['inference_data'], var_name="log_lik", pointwise=True)
res

/home/br1/PythonProjects/PythonStats/.venv/lib/python3.10/site-packages/arviz/stats/stats.py:803: UserWarning: Estimated shape parameter of Pareto distribution is greater than 0.7 for one or more samples. You should consider using a more robust model, this is because importance sampling is less likely to work well if the marginal posterior and LOO posterior are very different. This is more likely to happen with a non-robust model and highly influential observations.
  warnings.warn(

Computed from 10000 posterior samples and 17191 observations log-likelihood matrix.

         Estimate       SE
elpd_loo  7196.51   155.72
p_loo      260.37        -

There has been a warning during the calculation. Please check the results.
------

Pareto k diagnostic values:
                         Count   Pct.
(-Inf, 0.5]   (good)     17096   99.4%
 (0.5, 0.7]   (ok)          82    0.5%
   (0.7, 1]   (bad)         13    0.1%
   (1, Inf)   (very bad)     0    0.0%

aux_plt = az.plot_khat(res, show_bins=True, figsize=(7,5))

res = az.waic(ritorno['inference_data'], var_name="log_lik")
res

/home/br1/PythonProjects/PythonStats/.venv/lib/python3.10/site-packages/arviz/stats/stats.py:1645: UserWarning: For one or more samples the posterior variance of the log predictive densities exceeds 0.4. This could be indication of WAIC starting to fail. 
See http://arxiv.org/abs/1507.04544 for details
  warnings.warn(

Computed from 10000 posterior samples and 17191 observations log-likelihood matrix.

          Estimate       SE
elpd_waic  7208.48   155.58
p_waic      248.40        -

There has been a warning during the calculation. Please check the results.

plt_aux = az.plot_trace(ritorno['inference_data'], var_names=['phi','theta','hyper_gamma_phi_m','hyper_gamma_phi_s'], divergences=True, figsize=(15,15))

plt_aux = az.plot_trace(ritorno['inference_data'], var_names=['mh_sigma','sh_sigma','sigma'], divergences=True, figsize=(15,10))

plt_aux = az.plot_trace(ritorno['inference_data'], var_names=['y_0','betas','c'], divergences=True, figsize=(15,10))

df_risultati = az.summary(ritorno['inference_data'], var_names=['betas'], hdi_prob=0.95)
df_risultati.loc[:,'Name'] = ['Altitude','Type_Suburban','Type_Rural','Zoning_2','Zoning_3'] #????????? Come sono state fatte le covariate?
df_risultati.set_index('Name', inplace=True)
df_risultati

	mean	sd	hdi_2.5%	hdi_97.5%	mcse_mean	mcse_sd	ess_bulk	ess_tail	r_hat
Name
Altitude	-0.092	0.009	-0.108	-0.075	0.0	0.0	13429.0	8773.0	1.0
Type_Suburban	0.041	0.022	0.000	0.086	0.0	0.0	18357.0	7663.0	1.0
Type_Rural	0.084	0.027	0.031	0.136	0.0	0.0	18160.0	8105.0	1.0
Zoning_2	0.030	0.024	-0.018	0.077	0.0	0.0	12391.0	7872.0	1.0
Zoning_3	0.012	0.023	-0.033	0.058	0.0	0.0	12727.0	8524.0	1.0

import time
from matplotlib.patches import Rectangle

plt.figure(figsize=(10,3))
ax = plt.subplot(1,1,1)
plt.xlim(left = -0.3, right = 0.3)
pos_rows = np.arange(0,len(df_risultati.index))
plt.ylim(bottom=-0.75, top=pos_rows[-1]+1.0)
i = 1
for row in df_risultati.index:    
    lci = df_risultati.loc[row, 'hdi_2.5%']
    rci = df_risultati.loc[row, 'hdi_97.5%']
    len_ci = rci - lci
    if lci*rci <= 0:
        col = 'red'
    else:
        col = 'blue'
    ax.add_patch(Rectangle((lci, pos_rows[-i] - 0.75/2), len_ci, 0.75,
             edgecolor = 'grey',
             facecolor = col,
             fill=True,
             lw=1))
    i += 1
    
    
plt.axvline(0, 0, pos_rows[-1]+0.75, linestyle='--', color='black')
plt.grid(axis='both')
plt.yticks(ticks=np.flip(pos_rows), labels=df_risultati.index.to_list())
plt.xticks(ticks=np.linspace(-0.3,0.3,13))
plt.show()

aux_shape = ritorno['inference_data'].posterior.gamma_th.values.shape
sns.kdeplot(ritorno['inference_data'].posterior.gamma_th.values.reshape((aux_shape[0]*aux_shape[1])), legend=False)

<AxesSubplot: ylabel='Density'>

aux_shape = ritorno['inference_data'].posterior.phi.values.shape
sns.kdeplot(ritorno['inference_data'].posterior.phi.values.reshape((aux_shape[0]*aux_shape[1],aux_shape[2],aux_shape[3]))[:,0,:], legend=False)

<AxesSubplot: ylabel='Density'>

aux_shape = ritorno['inference_data'].posterior.phi.values.shape
sns.kdeplot(ritorno['inference_data'].posterior.phi.values.reshape((aux_shape[0]*aux_shape[1],aux_shape[2],aux_shape[3]))[:,1,:], legend=False)

<AxesSubplot: ylabel='Density'>

aux_shape = ritorno['inference_data'].posterior.c.values.shape
sns.kdeplot(ritorno['inference_data'].posterior.c.values.reshape((aux_shape[0]*aux_shape[1],aux_shape[2])), legend=False)

<AxesSubplot: ylabel='Density'>

aux_shape = ritorno['inference_data'].posterior.annual_mean.values.shape
sns.kdeplot(ritorno['inference_data'].posterior.annual_mean.values.reshape((aux_shape[0]*aux_shape[1],aux_shape[2])), legend=False)

<AxesSubplot: ylabel='Density'>

aux_shape = ritorno['inference_data'].posterior.annual_median.values.shape
sns.kdeplot(ritorno['inference_data'].posterior.annual_median.values.reshape((aux_shape[0]*aux_shape[1],aux_shape[2])), legend=False)

<AxesSubplot: ylabel='Density'>

aux_shape = ritorno['inference_data'].posterior.annual_max.values.shape
sns.kdeplot(ritorno['inference_data'].posterior.annual_max.values.reshape((aux_shape[0]*aux_shape[1],aux_shape[2])), legend=False)

<AxesSubplot: ylabel='Density'>

valori_31_dic = []
for dizionario in ritorno['missing']:
    if dizionario['Data'] == '2018-12-31':
        valori_31_dic.append(dizionario['Samples'])
        
res_plot = sns.kdeplot(valori_31_dic, legend=False)

aux_shape = ritorno['inference_data'].posterior.annual_days_over_threshold.values.shape
sns.kdeplot(ritorno['inference_data'].posterior.annual_days_over_threshold.values.reshape((aux_shape[0]*aux_shape[1],aux_shape[2])), legend=False)

<AxesSubplot: ylabel='Density'>

pd.set_option('display.max_columns', 500)
aux_results = pd.DataFrame(ritorno['inference_data'].posterior.annual_days_over_threshold.values.reshape((aux_shape[0]*aux_shape[1],aux_shape[2])), columns=df.columns).describe()
aux_results.sort_values('mean', axis=1, ascending=True)

	FEBBIO	CASTELLUCCIO	CORTE BRUGNATELLA	SAN LEO	BADIA	LUGAGNANO	VERUCCHIO	BESENZONE	SAN FRANCESCO	CASTELLARANO	GHERARDI	S. ROCCO	DELTA CERVIA	SAVIGNANO DI RIGO	VIA CHIARINI	GAVELLO	SAN PIETRO CAPOFIUME	PARCO MONTECUCCO	SARAGAT	S. LAZZARO	PARCO EDILCARANI	CITTADELLA	SAVIGNANO	GIARDINI MARGHERITA	REMESINA	MARECCHIA	FRANCHINI-ANGELONI	CENTO	PARCO BERTOZZI	VILLA FULVIA	CAORLE	PARCO RESISTENZA	PARCO FERRARI	SAN LAZZARO	PORTA SAN FELICE	MONTEBELLO	GIORDANI-FARNESE	MALCANTONE	TIMAVO	DE AMICIS	ROMA	PARADIGNA	GERBIDO	BOGOLESE	FLAMINIA	GIARDINI	CENO	ZALAMELLA	ISONZO
count	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.00000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.00000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.00000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000
mean	0.088900	0.242000	0.518000	1.658400	11.020100	14.089700	16.837300	18.176100	19.480300	20.329400	21.36450	21.736300	21.887300	22.496700	23.863300	23.898600	24.915300	25.972200	26.023300	27.877600	28.077500	28.635400	28.987000	29.153200	30.280900	30.60900	30.747000	30.889900	31.629800	32.483600	32.745600	34.145500	35.006800	36.476500	38.31910	39.277100	40.215400	40.230500	40.491100	43.092000	43.377800	45.209100	45.802600	46.167800	46.348100	46.421700	47.087500	47.689500	48.627400
std	0.354414	0.677999	1.108692	2.406136	7.924391	9.687045	11.346422	12.039817	11.469236	10.896815	14.15657	12.803099	12.822412	16.052176	14.165009	13.684862	15.842778	15.025372	14.112033	14.985682	13.541914	14.899163	15.550164	16.129642	15.542392	16.62208	16.944359	16.367558	18.132077	17.190151	18.233303	17.501079	17.200521	18.957716	18.74206	18.259092	18.930273	20.499302	19.149877	20.189678	20.412431	21.212887	22.656124	21.263044	21.166925	19.443422	22.314534	21.911634	20.719428
min	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.00000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	1.000000	0.000000	0.000000	1.000000	0.000000	0.00000	1.000000	0.000000	0.000000	0.000000	0.000000	1.000000	0.000000	1.000000	1.00000	1.000000	1.000000	1.000000	1.000000	0.000000	2.000000	0.000000	3.000000	1.000000	1.000000	4.000000	2.000000	2.000000	2.000000
25%	0.000000	0.000000	0.000000	0.000000	5.000000	7.000000	9.000000	10.000000	11.000000	12.000000	11.00000	13.000000	13.000000	11.000000	14.000000	14.000000	14.000000	15.000000	16.000000	17.000000	18.000000	18.000000	18.000000	17.000000	19.000000	19.00000	18.000000	19.000000	19.000000	20.000000	20.000000	21.000000	23.000000	23.000000	25.00000	26.000000	27.000000	25.000000	27.000000	28.000000	28.000000	30.000000	29.000000	31.000000	31.000000	32.000000	31.000000	32.000000	34.000000
50%	0.000000	0.000000	0.000000	1.000000	9.000000	12.000000	14.000000	16.000000	17.000000	18.000000	18.00000	19.000000	19.000000	19.000000	21.000000	22.000000	22.000000	23.000000	24.000000	25.000000	26.000000	26.000000	26.000000	26.000000	28.000000	28.00000	28.000000	28.000000	28.000000	29.500000	29.000000	32.000000	32.000000	34.000000	36.00000	37.000000	37.000000	37.000000	38.000000	40.000000	40.000000	42.000000	42.000000	43.000000	43.000000	44.000000	44.000000	45.000000	46.000000
75%	0.000000	0.000000	1.000000	2.000000	15.000000	19.000000	22.000000	24.000000	26.000000	26.000000	28.00000	28.000000	29.000000	30.000000	31.000000	31.000000	33.000000	33.000000	34.000000	36.000000	35.000000	37.000000	37.000000	38.000000	39.000000	39.00000	40.000000	40.000000	41.000000	42.000000	43.000000	44.000000	45.000000	47.000000	49.00000	49.000000	51.000000	52.000000	51.000000	54.000000	55.000000	57.000000	59.000000	58.000000	58.000000	58.000000	60.000000	60.000000	61.000000
max	5.000000	11.000000	15.000000	27.000000	77.000000	99.000000	115.000000	128.000000	101.000000	85.000000	119.00000	115.000000	100.000000	180.000000	110.000000	134.000000	156.000000	139.000000	123.000000	112.000000	113.000000	120.000000	139.000000	134.000000	145.000000	128.00000	148.000000	131.000000	157.000000	156.000000	152.000000	140.000000	155.000000	169.000000	157.00000	134.000000	156.000000	159.000000	157.000000	158.000000	204.000000	155.000000	173.000000	179.000000	173.000000	141.000000	174.000000	188.000000	147.000000

pd.set_option('display.max_columns', 500)
aux_results = pd.DataFrame(ritorno['inference_data'].posterior.is_over_daily_limit.values.reshape((aux_shape[0]*aux_shape[1],aux_shape[2])), columns=df.columns).describe()
aux_results.sort_values('mean', axis=1, ascending=True)

	CASTELLUCCIO	FEBBIO	CORTE BRUGNATELLA	SAN LEO	BADIA	LUGAGNANO	VERUCCHIO	BESENZONE	SAN FRANCESCO	CASTELLARANO	S. ROCCO	DELTA CERVIA	GHERARDI	GAVELLO	VIA CHIARINI	SAVIGNANO DI RIGO	SAN PIETRO CAPOFIUME	SARAGAT	PARCO MONTECUCCO	PARCO EDILCARANI	S. LAZZARO	CITTADELLA	SAVIGNANO	GIARDINI MARGHERITA	REMESINA	MARECCHIA	FRANCHINI-ANGELONI	CENTO	PARCO BERTOZZI	VILLA FULVIA	CAORLE	PARCO RESISTENZA	PARCO FERRARI	SAN LAZZARO	PORTA SAN FELICE	MONTEBELLO	MALCANTONE	GIORDANI-FARNESE	TIMAVO	DE AMICIS	ROMA	GERBIDO	PARADIGNA	BOGOLESE	CENO	FLAMINIA	ZALAMELLA	GIARDINI	ISONZO
count	10000.0	10000.0	10000.0	10000.0	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.00000	10000.00000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.00000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.00000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000
mean	0.0	0.0	0.0	0.0	0.012700	0.034300	0.065300	0.088000	0.088500	0.09430	0.12990	0.136100	0.140200	0.169500	0.176500	0.177800	0.210500	0.216100	0.217500	0.248000	0.258600	0.276900	0.288000	0.291100	0.31400	0.316900	0.324600	0.328900	0.341000	0.363300	0.366400	0.407700	0.427600	0.45310	0.501000	0.522000	0.532300	0.540700	0.549200	0.598500	0.601900	0.627300	0.633700	0.653200	0.660500	0.660900	0.678800	0.679400	0.715700
std	0.0	0.0	0.0	0.0	0.111982	0.182008	0.247067	0.283309	0.284035	0.29226	0.33621	0.342912	0.347212	0.375212	0.381264	0.382363	0.407684	0.411604	0.412566	0.431873	0.437887	0.447489	0.452854	0.454292	0.46414	0.465291	0.468249	0.469837	0.474069	0.480974	0.481845	0.491431	0.494755	0.49782	0.500024	0.499541	0.498981	0.498366	0.497598	0.490226	0.489531	0.483547	0.481817	0.475975	0.473563	0.473428	0.466961	0.466731	0.451103
min	0.0	0.0	0.0	0.0	0.000000	0.000000	0.000000	0.000000	0.000000	0.00000	0.00000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.00000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.00000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000
25%	0.0	0.0	0.0	0.0	0.000000	0.000000	0.000000	0.000000	0.000000	0.00000	0.00000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.00000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.00000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000
50%	0.0	0.0	0.0	0.0	0.000000	0.000000	0.000000	0.000000	0.000000	0.00000	0.00000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.00000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.00000	1.000000	1.000000	1.000000	1.000000	1.000000	1.000000	1.000000	1.000000	1.000000	1.000000	1.000000	1.000000	1.000000	1.000000	1.000000
75%	0.0	0.0	0.0	0.0	0.000000	0.000000	0.000000	0.000000	0.000000	0.00000	0.00000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	1.000000	1.000000	1.000000	1.000000	1.00000	1.000000	1.000000	1.000000	1.000000	1.000000	1.000000	1.000000	1.000000	1.00000	1.000000	1.000000	1.000000	1.000000	1.000000	1.000000	1.000000	1.000000	1.000000	1.000000	1.000000	1.000000	1.000000	1.000000	1.000000
max	0.0	0.0	0.0	0.0	1.000000	1.000000	1.000000	1.000000	1.000000	1.00000	1.00000	1.000000	1.000000	1.000000	1.000000	1.000000	1.000000	1.000000	1.000000	1.000000	1.000000	1.000000	1.000000	1.000000	1.00000	1.000000	1.000000	1.000000	1.000000	1.000000	1.000000	1.000000	1.000000	1.00000	1.000000	1.000000	1.000000	1.000000	1.000000	1.000000	1.000000	1.000000	1.000000	1.000000	1.000000	1.000000	1.000000	1.000000	1.000000

aux_results.loc['mean',:].plot(kind='barh', figsize=(15,15), xticks=np.linspace(0,1,11))
#res = plt.xticks(rotation = 30)

<AxesSubplot: >

pd.set_option('display.max_columns', 500)
aux_results = pd.DataFrame(ritorno['inference_data'].posterior.is_over_annual_limit.values.reshape((aux_shape[0]*aux_shape[1],aux_shape[2])), columns=df.columns).describe()
aux_results.sort_values('mean', axis=1, ascending=True)

	CASTELLUCCIO	GHERARDI	SARAGAT	S. LAZZARO	S. ROCCO	DELTA CERVIA	PARCO EDILCARANI	CASTELLARANO	LUGAGNANO	VIA CHIARINI	CITTADELLA	CENTO	SAN FRANCESCO	FEBBIO	CORTE BRUGNATELLA	SAN LEO	VERUCCHIO	MARECCHIA	BADIA	BESENZONE	PARCO MONTECUCCO	REMESINA	SAVIGNANO	GAVELLO	PARCO FERRARI	GIARDINI	GIARDINI MARGHERITA	PARCO RESISTENZA	MONTEBELLO	FRANCHINI-ANGELONI	SAN PIETRO CAPOFIUME	SAVIGNANO DI RIGO	CAORLE	VILLA FULVIA	GIORDANI-FARNESE	PARCO BERTOZZI	TIMAVO	PORTA SAN FELICE	SAN LAZZARO	ISONZO	DE AMICIS	ROMA	MALCANTONE	FLAMINIA	PARADIGNA	ZALAMELLA	BOGOLESE	CENO	GERBIDO
count	10000.0	10000.0	10000.0	10000.0	10000.0	10000.0	10000.0	10000.0	10000.0	10000.0	10000.0	10000.0	10000.0	10000.0	10000.0	10000.0	10000.0	10000.0	10000.0	10000.0000	10000.0000	10000.0000	10000.0000	10000.0000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.00000	10000.00000	10000.000000	10000.00000	10000.000000	10000.000000	10000.000000	10000.000000	10000.000000	10000.00000	10000.000000	10000.000000
mean	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0001	0.0001	0.0001	0.0001	0.0001	0.000200	0.000300	0.000300	0.000300	0.000300	0.000300	0.000300	0.000400	0.000500	0.000900	0.000900	0.001000	0.001000	0.00110	0.00110	0.001200	0.00160	0.001700	0.002300	0.002600	0.003000	0.003700	0.00380	0.004600	0.005100
std	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0100	0.0100	0.0100	0.0100	0.0100	0.014141	0.017319	0.017319	0.017319	0.017319	0.017319	0.017319	0.019997	0.022356	0.029988	0.029988	0.031609	0.031609	0.03315	0.03315	0.034622	0.03997	0.041198	0.047906	0.050926	0.054693	0.060718	0.06153	0.067671	0.071236
min	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0000	0.0000	0.0000	0.0000	0.0000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.00000	0.00000	0.000000	0.00000	0.000000	0.000000	0.000000	0.000000	0.000000	0.00000	0.000000	0.000000
25%	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0000	0.0000	0.0000	0.0000	0.0000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.00000	0.00000	0.000000	0.00000	0.000000	0.000000	0.000000	0.000000	0.000000	0.00000	0.000000	0.000000
50%	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0000	0.0000	0.0000	0.0000	0.0000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.00000	0.00000	0.000000	0.00000	0.000000	0.000000	0.000000	0.000000	0.000000	0.00000	0.000000	0.000000
75%	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0000	0.0000	0.0000	0.0000	0.0000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.00000	0.00000	0.000000	0.00000	0.000000	0.000000	0.000000	0.000000	0.000000	0.00000	0.000000	0.000000
max	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	1.0000	1.0000	1.0000	1.0000	1.0000	1.000000	1.000000	1.000000	1.000000	1.000000	1.000000	1.000000	1.000000	1.000000	1.000000	1.000000	1.000000	1.000000	1.00000	1.00000	1.000000	1.00000	1.000000	1.000000	1.000000	1.000000	1.000000	1.00000	1.000000	1.000000

aux_results.loc['mean',:].plot(kind='barh', figsize=(15,15), xticks=np.linspace(0,1,11))
#res = plt.xticks(rotation = 30)

<AxesSubplot: >

y_post_pred = ritorno['inference_data'].posterior.y_post_pred.to_numpy()
y_post_pred = y_post_pred.reshape(y_post_pred.shape[0]*y_post_pred.shape[1], y_post_pred.shape[2], y_post_pred.shape[3])
y_post_pred.shape

(10000, 365, 49)

import ipywidgets as widgets
from ipywidgets import HBox, VBox
from IPython.display import display

@widgets.interact(stazione = df.columns)
def f(stazione):
    col_map = sns.light_palette((20, 75, 70), input='husl', as_cmap=True)

    plt.figure(figsize=(14,8))
    ax = plt.subplot(1,1,1)

    idx_stazione = df.columns.to_list().index(stazione)
    
    """
    sns.lineplot(np.transpose(y_post_pred[0:50,:,idx_stazione]))
    """

    lower_lim = np.zeros(len(df.index))
    upper_lim = np.zeros(len(df.index))    
    for i in range(len(df.index)):
        lower_lim[i] = np.percentile(y_post_pred[:,i,idx_stazione],2.5)
        upper_lim[i] = np.percentile(y_post_pred[:,i,idx_stazione],97.5)
    
    
#    sns.lineplot(lower_lim)
#    sns.lineplot(upper_lim)
    
    for i in range(len(df.index)):
        ax.add_patch(mplt.patches.Rectangle((i,lower_lim[i]),1,upper_lim[i]-lower_lim[i], fill=True, color=col_map(180)))

    sns.lineplot(df[stazione])    
    sns.lineplot(np.mean(y_post_pred[:,:,idx_stazione], axis=0))


    index = 0
    for line in ax.get_lines():
        if index == 0:
            col_map = sns.dark_palette((230,90,65), input='husl', as_cmap=True)
            line.set_c(col_map(155))
        else:
            col_map = sns.dark_palette((10,90,65), input='husl', as_cmap=True)
            line.set_c(col_map(175))
        index += 1
        
    plt.ylim(bottom=0,top=2.5)

    primo_giorno = datetime.date(2018,1,1)
    date_da_segnare = []
    date_da_segnare_posizioni = []

    for i in range(12):
        date_da_segnare.append(datetime.date(2018,i+1,1))
        date_da_segnare_posizioni.append((date_da_segnare[2*i] - primo_giorno).days)
        date_da_segnare.append(datetime.date(2018,i+1,15))
        date_da_segnare_posizioni.append((date_da_segnare[2*i+1] - primo_giorno).days)
        date_da_segnare[2*i] = date_da_segnare[2*i].isoformat()
        date_da_segnare[2*i+1] = date_da_segnare[2*i+1].isoformat()

    plt.xticks(date_da_segnare_posizioni,date_da_segnare,rotation=65)
    plt.tick_params(
        axis='x',          # changes apply to the x-axis
        which='both',      # both major and minor ticks are affected
        bottom=True,      # ticks along the bottom edge are off
        top=False,         # ticks along the top edge are off
        labelbottom=True)
    plt.grid()
    
    """
    ax.get_legend().remove()
    col_vals = np.linspace(1,255,num=len(df.columns))
    index = 0
    for line in ax.get_lines():
        if(index == len(ax.get_lines()) - 1):
            col_map = sns.dark_palette((230,90,65), input='husl', as_cmap=True)
            line.set_c(col_map(175))
            continue
        if(index == len(ax.get_lines()) - 2):
            col_map = sns.dark_palette((120,90,65), input='husl', as_cmap=True)
            line.set_c(col_map(175))
            index += 1
            continue

        #line.set_c(col_map(int(np.round(col_vals[index]))))
        line.set_c(col_map(220))
        line.set_alpha(0.3)
        index += 1
    """
    
    

    plt.show()

interactive(children=(Dropdown(description='stazione', options=('CASTELLUCCIO', 'FEBBIO', 'CORTE BRUGNATELLA',…

def f(stazione):
    col_map = sns.light_palette((20, 75, 70), input='husl', as_cmap=True)

    plt.figure(figsize=(14,8))
    ax = plt.subplot(1,1,1)

    idx_stazione = df.columns.to_list().index(stazione)
    
    """
    sns.lineplot(np.transpose(y_post_pred[0:50,:,idx_stazione]))
    """

    lower_lim = np.zeros(len(df.index))
    upper_lim = np.zeros(len(df.index))    
    for i in range(len(df.index)):
        lower_lim[i] = np.percentile(y_post_pred[:,i,idx_stazione],2.5)
        upper_lim[i] = np.percentile(y_post_pred[:,i,idx_stazione],97.5)
    
    
#    sns.lineplot(lower_lim)
#    sns.lineplot(upper_lim)
    
    for i in range(len(df.index)):
        ax.add_patch(mplt.patches.Rectangle((i,lower_lim[i]),1,upper_lim[i]-lower_lim[i], fill=True, color=col_map(180)))

    sns.lineplot(df[stazione])    
    sns.lineplot(np.mean(y_post_pred[:,:,idx_stazione], axis=0))


    index = 0
    for line in ax.get_lines():
        if index == 0:
            col_map = sns.dark_palette((230,90,65), input='husl', as_cmap=True)
            line.set_c(col_map(155))
        else:
            col_map = sns.dark_palette((10,90,65), input='husl', as_cmap=True)
            line.set_c(col_map(175))
        index += 1
        
    plt.ylim(bottom=0,top=2.5)

    primo_giorno = datetime.date(2018,1,1)
    date_da_segnare = []
    date_da_segnare_posizioni = []

    for i in range(12):
        date_da_segnare.append(datetime.date(2018,i+1,1))
        date_da_segnare_posizioni.append((date_da_segnare[2*i] - primo_giorno).days)
        date_da_segnare.append(datetime.date(2018,i+1,15))
        date_da_segnare_posizioni.append((date_da_segnare[2*i+1] - primo_giorno).days)
        date_da_segnare[2*i] = date_da_segnare[2*i].isoformat()
        date_da_segnare[2*i+1] = date_da_segnare[2*i+1].isoformat()

    plt.xticks(date_da_segnare_posizioni,date_da_segnare,rotation=65)
    plt.tick_params(
        axis='x',          # changes apply to the x-axis
        which='both',      # both major and minor ticks are affected
        bottom=True,      # ticks along the bottom edge are off
        top=False,         # ticks along the top edge are off
        labelbottom=True)
    plt.grid()
    
    """
    ax.get_legend().remove()
    col_vals = np.linspace(1,255,num=len(df.columns))
    index = 0
    for line in ax.get_lines():
        if(index == len(ax.get_lines()) - 1):
            col_map = sns.dark_palette((230,90,65), input='husl', as_cmap=True)
            line.set_c(col_map(175))
            continue
        if(index == len(ax.get_lines()) - 2):
            col_map = sns.dark_palette((120,90,65), input='husl', as_cmap=True)
            line.set_c(col_map(175))
            index += 1
            continue

        #line.set_c(col_map(int(np.round(col_vals[index]))))
        line.set_c(col_map(220))
        line.set_alpha(0.3)
        index += 1
    """
    
    

    plt.show()

for stazione in df.columns:
    f(stazione)